A theoretical modeling framework for motile and colonial harmful algae

Abstract Climate change is leading to an increase in severity, frequency, and distribution of harmful algal blooms across the globe. For many harmful algae species in eutrophic lakes, the formation of such blooms is controlled by three factors: the lake hydrodynamics, the vertical motility of the algae organisms, and the ability of the organisms to form colonies. Here, using the common cyanobacterium Microcystis aeruginosa as an example, we develop a model that accounts for both vertical transport and colony dynamics. At the core of this treatment is a model for aggregation. For this, we used Smoluchowski dynamics containing parameters related to Brownian motion, turbulent shear, differential setting, and cell‐to‐cell adhesion. To arrive at a complete description of bloom formation, we place the Smoluchowski treatment as a reaction term in a set of one‐dimensional advection‐diffusion equations, which account for the vertical motion of the algal cells through molecular and turbulent diffusion and self‐regulating buoyant motion. Results indicate that Smoluchowski aggregation qualitatively describes the colony dynamics of M. aeruginosa. Further, the model demonstrates wind‐induced mixing is the dominant aggregation process, and the rate of aggregation is inversely proportional to algal concentration. Because blooms of Microcystis typically consist of large colonies, both of these findings have direct consequences to harmful algal bloom formation. While the theoretical framework outlined in this manuscript was derived for M. aeruginosa, both motility and colony formation are common among bloom‐forming algae. As such, this coupling of vertical transport and colony dynamics is a useful step for improving forecasts of surface harmful algal blooms.


| INTRODUC TI ON
Microcystis aeruginosa is a common toxin-producing cyanobacterium capable of forming harmful algal blooms (HABs). HABs threaten both ecological and public health, and they are expected to increase in distribution, frequency, and severity as a result of climate change (O'neil et al., 2012). Predicting the timing of bloom formation has been challenging, but researchers in the field have reached consensus on general trends leading up to a HAB. A study of the recordbreaking Lake Erie algae bloom of 2011 determined that-in addition to excessive nutrient loading-quiescent meteorological conditions allowed the bloom to form and proliferate to such a massive extent (Michalak et al., 2013), a finding that has been corroborated in many subsequent studies of cyanobacteria HABs (Wells et al., 2015). Using a Bayesian biophysical model with a high-frequency dataset, Del Giudice et al. (2021) were able to quantitatively show that quiescent conditions are not enough: High surface water temperatures and high irradiation are also necessary for bloom formation. Recently, it has been suggested that vertical heterogeneity (i.e., subsurface peaks) of M. aeruginosa concentration is an important precursor to Microcystis surface bloom formation (Seegers et al., 2015;Xiao et al., 2018;Wilkinson et al., 2019;Taylor et al., 2021). Therefore, it is reasonable to assume improving models for the drivers of M. aeruginosa vertical distributions will likely lead to improved predictions of HAB timing.
There are two key traits related to the ubiquity of M. aeruginosa: vertical motility and colony formation. Vertical motility is achieved through algal cell buoyancy regulation via intracellular gas vesicles.
Under low levels of mixing, M. aeruginosa sinks to lower light intensities during the day and floats towards the water surface at night, although a critical water temperature threshold must be reached in order for cells to regain buoyancy (Ibelings et al., 1991;Thomas & Walsby, 1985, 1986. Once that threshold is reached, increasing temperature increases buoyant velocity (You et al., 2018). Vertical motility gives M. aeruginosa a particular advantage in stratified lake environments. Stratified lakes are characterized by three distinct layers: The epilimnion or surface mixed layer is the hot, well-mixed surface layer; the hypolimnion is the cold, well-mixed bottom layer; and the metalimnion is the intermediate layer of steep temperature gradient connecting the epilimnion to the hypolimnion. Using the three-dimensional ecological-hydrodynamic modeling software ELCOM-CAEDYM, Chung et al. (2014) were able to demonstrate a shallow mixed layer depth (close to the photic depth) favored buoyant cyanobacteria dominance, indicating lake thermal structure controls algal population dynamics.
Colony dynamics remain rather illusive, but colonies have been demonstrated to form in the presence of grazers, low to medium turbulence, and low nutrient conditions. Colonies formed by reproduction and growth tend to be compact, whereas colonies that form by collisions tend to be fractal. There is also a well-documented progression from a unicellular morphology in the spring to a fractal colonial morphology in the summer (Xiao et al., 2018). In a field study, Cao and Yang (2010) found that large colonies (greater than 20 cells per colony) did not appear until May but composed 90% of cells in a June surface bloom. They also calculated the mean number of cells in the surface bloom to be about 120 cells/colony. Between field work and experiments, Qin et al. (2018) found that wind promotes aggregation, creating heterogeneous size distributions in Microcystis populations.
There are two threads of previous models to follow. There are models that describe aggregation processes of phytoplankton, and there are models that describe the vertical motility of M. aeruginosa. To describe the aggregation processes of phytoplankton, models use Smoluchowski aggregation terms (Ackleh & Miller, 2018;Jackson, 1990;Smoluchowski, 1917). Because these models typically have applications in wastewater treatment or marine snow, the only transport considered is the loss of aggregates via sinking out of the surface mixed layer (Engel et al., 2004;Lee et al., 2000;Teh et al., 2016).
Early models of Microcystis motility use light intensity as a driver of changes in individual cell density-high light intensities lead to an increase in cell density, whereas low light intensities lead to a decrease in cell density. The buoyant velocity of cells is then calculated through a modified Stokes settling velocity that is governed by the difference between algal cell density and the surrounding water density (Wallace et al., 2000). Turbulent transport has since been incorporated into these models (Medrano et al., 2013;Zhu et al., 2018). By combining their model with principal component analysis, Feng et al. (2018) demonstrated that turbulence-induced mixing explained over half of the variability of early surface bloom formation, and that buoyancy regulation was more important for bloom maintenance and formation of late-season blooms. Although the transport of different (fixed) colony sizes is investigated in the aforementioned Microcystis motility models, they do not incorporate aggregation dynamics, despite the well-documented progression from unicellular to colonial morphologies.
In a previous field study, statistical methods were used to elucidate the reliance of Microcystis-dominated algal vertical distributions on lake thermal stratification variables (Taylor et al., 2021).
Following the protocol discussed in Vinatier et al. (2011), which suggests using statistical and mechanistic models in an iterative manner to uncover forcings of spatial heterogeneity, we propose a mechanistic model to analyze the effects of hydrodynamic and biological processes underlying the spatial patterns observed in the previous field study. The primary objective of this model is not to replicate exact field observations but to instead generate hypotheses for the biophysical drivers of general field trends and observations. To this end, we couple algal cell aggregation dynamics with algal motility in a system of one-dimensional partial differential equations that capture lake hydrodynamics to investigate the role of the colony and motility dynamics on M. aeruginosa surface bloom formation.

| Aggregation preliminaries
In the absence of any advective or diffusive transport, discrete aggregation dynamics can be described by the Smoluchowski coagulation model (Smoluchowski, 1917): where n k (z,t) is the concentration of an aggregate of size k, α(i,j) is the sticking probability and β(i,j) is referred to as the aggregation, or coagulation, kernel of particles of size i and j (Figure 1). Occasionally the product of α(i,j) and β(i,j) is referred to as the aggregation kernel, instead of just β(i,j). We leave the two parameters decoupled mainly for the sake of visualizing the process (Figure 1) but also to conceptually differentiate the hydrodynamic drivers of β(i,j) (Equations 2-5) from the biological drivers of α(i,j) (Section 2.2.2). The first term on the right-hand side describes the formation of a k-sized aggregation, whereas the second term on the right-hand side describes the loss of a k-sized aggregation through the formation of a k + i-sized aggregate. An infinitely-sized particle represents a loss of mass due to gelation. Equation 1 has had far-reaching applications in addition to phytoplankton modeling, from aerosols to random graph theory and polymerization to planet formation (Aldous, 1999).
While analytical solutions exist for some simple aggregation kernels ( (i, j) ∼ 1, (i, j) ∼ i + j, and (i, j) ∼ ij), realistic aggregation kernels are rarely analytically tractable. In the present context, β(i,j) is calculated as the sum of aggregation kernels for Brownian motion, β Br (i,j,z), turbulent shear, β TS (i,j,z), and differential settling, β DS (i,j,z), each, respectively, defined as (Ackleh & Miller, 2018;Thomas et al., 1999) and such that where T(z) is the water temperature (K), k B is Boltzmann's constant (1.38 × 10 −23 m 2 kg s −2 K −1 ), μ(z) is the dynamic viscosity of water (kg/ m/s), G(z) = 1 2 is the turbulent shear rate (1/s), ϵ(z) is the rate of turbulent kinetic energy dissipation (m 2 /s 3 ), and (z) is the kinematic viscosity of water (m 2 /s). The equivalent spherical diameter of a colony of size i, d i (m), is given by where D f = 2.5 is the fractal dimension (Nakamura et al., 1993), d 0 = 5 μm is the diameter of a single cell of M. aeruginosa (Xiao et al., 2018), and ϕ is the colony porosity that linearly decreases from ϕ = 1 for single cells and ϕ = 0.2 for colonies of size k max (Medrano et al., 2013). Equation (2) is derived from thermodynamic principles of Brownian motion, Equation (3) defines the rate of collisions for sub-Kolmogorov particles in turbulent flow (i.e., the largest aggregate diameter is smaller than the length scale of the smallest turbulent eddies), and Equation (4) describes collisions as a result of differentsized aggregates moving at different velocities. Aggregation due to Brownian motion is typically much slower than aggregation due to turbulent shear, and aggregation due to differential settling will be large for aggregates of drastically different sizes but will be small for aggregates of close to the same size.
There are several assumptions of this formulation that should be addressed before continuing.
1. We assume diffusion-limited aggregation rather than reactionlimited aggregation, meaning the aggregation process will be limited by diffusion due to Brownian motion and not by the sticking probability of collisions. This is reasonable for colonyforming species of algae in a system where the domain size is much larger than the aggregate sizes.
2. We assume a maximum colony size, below which there will be no disaggregation-colonies cannot split up once formed. Effectively, we assume any colonies above the maximum colony size instanta- 3. We assume aggregates grow in size through particle collisions only. When aggregates consist of living organisms, it is possible for aggregates to increase in size through cell growth and reproduction in addition to particle collisions. However, it is hypothesized that the fractal colonies of M. aeruginosa are formed primarily through collisions, so we neglect aggregation due to cell growth (Xiao et al., 2018). 4. We assume aggregation is uniform over any given horizontal cross-section in order to facilitate the construction of a onedimensional model.

| The mathematical model
To Since M. aeruginosa buoyancy is largely mediated by light intensity, we must also construct diurnal light profiles. We generated surface light intensities, I 0 (t), by where I max is the maximum surface light intensity and D L is the photoperiod. To best replicate previous models, values of I max = 800 W/m 2 and D L = 16 h were chosen (Medrano et al., 2013 (Cao & Yang, 2010). Using this, we define a function that gives the Smoothed field data. Low wind profiles for (a) temperature, (b) turbulent dispersion coefficient, D Z , and (c) rate of turbulent kinetic energy dissipation, ϵ. high wind profiles for (d) temperature, (e) turbulent dispersion coefficient, D Z , and (f) rate of turbulent kinetic energy dissipation, ϵ. Note the differences in orders of magnitude for D Z and ϵ under low wind and high wind conditions. Low wind conditions roughly correspond to wind speeds of 2 m/s, whereas high wind conditions roughly correspond to wind speeds of 8 m/s New figure to elaborate on the sticking probability function. Sticking probability, α k , vs colony diameter, d k (μm), and colony size, k (cells/colony), where α k is defined by . Single cells will aggregate upon colliding 50% of the time, whereas colonies of size k = 95 cells/colony will always aggregate after collisions. Note that (i, j) = max i , j sticking probability of a colony of size k, α k = f(d k ), which achieves a minimum value of α k = 0.5 at d 1 = 5 μm and a maximum value of α k = 1 at d 95 = 125 μm ( Figure 3). To calculate the sticking probability for a collision between a colony of size i and size j, we define (i, j) = max i , j . Larger colonies will therefore be 'stickier' than small colonies, so more of their collisions will result in aggregation.
The buoyant velocity, w k , is calculated using subroutines described in previous models, which (i) relate light intensity to individual cell density, then (ii) relate individual cell density to colony density using the fractal dimension of M. aeruginosa aggregates, then (iii) use the colony density to calculate a modified Stoke's velocity (Medrano et al., 2013;Nakamura et al., 1993;Wallace et al., 2000) by where ρ k is the density of a colony of size k. Subroutine details to calculate ρ k can be found in Appendix A. We use the same equations and parameter values used in the work of Medrano et al. (2013), with a modification for the ratio of cell volume to colony volume that accounts for the fractal geometry of aggregates and the relationship between EPS content and colony size. We expect sinking during the day (positive w k ) and floating at night (negative w k ), although velocity magnitudes and general transport dynamics will vary across colony  Table 1 shows numerical parameter values used for all simulations. The time step, Δt, was chosen to be small enough to ensure the stability of the numerical scheme, and the grid cell width, Δz, was chosen to be small enough to minimize numerical dispersion of the upwind scheme while also maintaining stability. To address numerical dispersion, we tested the time to large colony appearance for the parameters described in Table 1 against a finer grid size. In the base case simulation, large colonies appear in 13.4 days; if we instead use Δz = 0.1 m (and a correspondingly smaller time step of Δt = 5 s), large colonies appear in 16.1 days. This three-day slowdown indicates that our scheme is not completely devoid of numerical dispersion. However, the goal of this manuscript is first and foremost to investigate the applicability of Smoluchowski aggregation to describe M. aeruginosa colony dynamics-not to solve the inverse problem of parameter estimation or make predictions with a real dataset. In this sense, we feel that our choices of space and time steps efficiently capture the correct physical behaviors and provide an appropriate order of magnitude prediction for the timing and appearance of large colony sizes.

| Appearance and distribution of colonies
We will start with the simplest simulation that still allows for the investigation of important model features: six weeks of a repeating photoperiod and constant lake thermal and hydrodynamic profiles ( Table 2). The repeating photoperiod is generated by Equations (11) (13)  Figure 2d-f. For the base case simulation, the lake thermal and hydrodynamic profiles represent high wind conditions. Field data indicate Microcystis can transition from a predominantly unicellular morphology to a predominantly colonial morphology over a monthly period (Cao & Yang, 2010;Xiao et al., 2018), so a six-week simulation time was chosen to ensure aggregation would be evident.
The model demonstrates small colonies will diffuse throughout the mixed layer (Figure 5a-c), but large colonies exhibit diurnal migrations to a depth with preferred low light intensity (Figure 5d,e).
In general, small colonies will lose mass as they aggregate into larger colonies, which gain mass. Medium-sized colonies never achieve high mass (Figure 5c,d), and colonies of size k = 101 appear before colonies of size k = 67. This indicates large colonies aggregate with each other faster than they aggregate with small colonies, a finding consistent with coagulation kinetic theory (Smit et al., 1994). The overall concentration profile, C(z,t) (Equation 10), is mostly influenced by large colonies by approximately the fifth week of simulation (Figure 5f).

| Factors affecting vertical distribution
While advection is negligible for single cells and small colonies, motility plays a key role in the vertical distribution of large-sized colonies ( Figure 6). The time it takes for large colonies to appear is approximately equivalent to whether advection is on or off, but the inclusion of motility allows the large colonies to migrate to a preferred depth of low light intensity (Figure 6a).
We also see changes in vertical distributions when we change wind conditions (Figure 7). During high wind conditions, small colonies become uniformly distributed throughout the mixed layer.  Figure 5e). In addition, wind also seems to significantly control the time it takes for colonies to appear. Synthesizing these results, high wind conditions lead to more medium-sized colonies, but they will be well-mixed throughout the surface mixed layer. On the other hand, low wind conditions lead to far fewer medium-sized colonies, but the colonies will be able to concentrate around a depth of preferred low light intensity.

| Factors affecting aggregation
There are few situations less likely to occur than 6 weeks of the exact same meteorological conditions on repeat, so we must explore how the model behaves under different conditions. To this end, let us define to be the total number of cells in a colony of size k. Since n k is a continuous variable and n k Δz is not necessarily greater than one, it is possible for N k < k. We are more concerned when colonies of various sizes appear at some comparative concentration value rather than the actual concentration, so N k (t) acts as a suitable marker for the appearance of colonies.
We can now rerun the simulation described in the previous Section 3.1 while changing one condition at a time to see how each individual change affects N k (t) for various colony sizes (Figures 8 and 9). Along with wind conditions, the speed of aggregation is highly sensitive to initial algal concentrations ( Figure 10). Let us define τ k to be the time such that N k (τ k ) = 1. As long as initial algal concentrations are greater than 1 × 10 7 cells/m 3 , then τ k is approximately inversely proportional to initial concentrations within the mixed layer, n 0 3. Aggregation is negligible during low wind conditions. 4. Intermittent wind conditions, which oscillate between high and low winds at some given frequency such that high wind conditions are achieved 50% of the time, slow the appearance of large colony sizes by a factor of two.
5. Above an initial algal concentration of 10 7 cells/m 3 , there is a power-law dependence between the time to appearance of large colonies and initial algal concentration. Figure 2d Cell counts, N k , were calculated by

F I G U R E 5 Concentration profiles over six weeks of simulation during high wind conditions (shown in
size to appear during low wind conditions was k = 3 cells/colony ( Figure 8a,b). Cutting the large wind events in half-either daily or hourly-slowed the appearance of the large-sized colonies by a factor of two (Figure 8c,d). This implies that the speed of aggregation is directly proportional to the duration of large wind events, causing relatively short-lived wind events to lead to rapid aggregation (recall the dependence of β(i,j,z) on the turbulent shear rate in Equation (3)).
When our model indicates aggregation is negligible for low wind conditions, it does not mean aggregation is not occurring.
Instead, this indicates that processes like light-driven motility are considerably more significant than aggregation during low wind conditions. While this may seem to disagree with the conclusions of Qin et al. (2018), which stated that low to medium turbulence is necessary to promote colony formation, this finding is actually just placing their experimental results in the context of a deep, dimictic lake. Small to moderate amounts of turbulence will in fact promote aggregation, but it does not do so at a rate that will lead to large colonies appearing in a six-week time frame. Furthermore, even in shallow Lake Taihu, Qin et al. (2018) measured a significant increase in average colony size over a short period of several days during a typhoon event with consistently high wind speeds, a result consistent with our findings.
This observation has profound consequences on the subsequent formation of surface blooms. Shortly after large wind events, the newly large colonies will be able to overcome turbulent mixing that the previously small colonies could not, leading to drastically different vertical transport results. Since blooms typically consist of large colonies (Cao & Yang, 2010;Wu et al., 2020;Zhu et al., 2014), this also means short periods of mixing via large wind events could act as a necessary precursor to surface harmful algal bloom formation as long as algal concentrations are high enough (see discussion below and Section 4.3 for more details). In a laboratory mesocosm experiment, Wu et al. (2019) found that increasing wind (up to 3.6 m/s) increased the volumetric median colony diameter at the water surface. Field experiments by Yang et al. (2020) found that intermittent wind-induced disturbance favored (i) larger colony sizes, (ii) higher biomass, and (iii) stronger dominance of Microcystis over constant quiescent or constant wind conditions. We believe this result agrees nicely with our conclusion that wind is necessary to promote aggregation, quiescent conditions are necessary for algal growth, and the combination of the two in subsequent order is a recipe for a harmful algal bloom.
In regards to the sensitivity of aggregation to the initial algal concentration, the inversely proportional relationship between algal concentration and time to large colony appearance, τ 101 , has been documented in previous studies of marine snow. Jackson (1990) found their large-sized colonies appeared within half a day of algal concentrations reaching 10 8 cells/m 3 , a rate in line with the results described in this manuscript ( Figure 10). We relate τ 101 to initial concentrations only, but that is simply because we have a conserved number of total cells in our system. If instead we had growth and/ or decay terms, we could track τ 101 as a function of instantaneous algal concentration. By maintaining conservation of mass, however, we can clearly see that any location in the water column with algal concentrations on the order of 10 7 cells/m 3 will take over 10 days to form large colonies, whereas locations with concentrations on the order of 10 8 cells/m 3 will have large colonies within a day.
Since higher densities would lead to increased collisions, this finding is unsurprising from a physical standpoint; however, it does provide some important biological modeling insight. Regardless of wind conditions, aggregation will be negligible until algal F I G U R E 9 Cell count, N k (t), of various colony sizes for (a) the base case simulation from Figure 5  ). With a starting concentration of 1 × 10 7 cells/m 3 , colonies of size k = 101 never appear within the 42-day simulation period concentration exceeds 10 7 cells/m 3 . After this threshold is reached, the rate of aggregation will increase as concentration increases. A large wind event later in the season-when algal concentrations are high-will therefore have dramatically different aggregation consequences than a large wind event in the beginning of the season, when algal concentrations are low. Further, nonuniform algal concentration profiles will lead to nonuniform aggregation. Any depth where there is a peak in algal concentration will also act as a hot spot for aggregation, leading to nonuniform colony size distributions within the water column.

| An evaluation of model assumptions
Before addressing the implications of these findings on harmful algal blooms, we must discuss how model assumptions may impact results. Let us start with our neglect of disaggregation and our limitation on maximum colony size. Large colonies (d k > 420 μm) would almost surely fragment under our high wind/strong turbulence conditions (O'Brien et al., 2004). The fact that turbulence also promotes aggregation through enhanced mixing represents a colony size trade-off. Turbulence causes the colony size distribution to skew towards the largest stable colony size, but the largest stable colony size decreases with increasing turbulence. If we were to allow larger colony sizes in the model, we would have to include fragmentation, a conclusion arrived at by Ackleh and Miller (2018)  Another constraint of this model is the restriction of algal growth, which is negligible over short timescales but significant over seasonal timescales. Recall this decision was made because M. aeruginosa colonies tend to be fractal in shape, and fractal aggregates are often the result of aggregation due to collisions instead of cell growth (Xiao et al., 2018). However, in experiments, Duan et al. (2018) found that Microcystis colony size significantly increased with increasing temperature. Although the aggregation kernel related to Brownian motion scales linearly with temperature (Equation (2)), this thermodynamic dependency alone cannot explain this variability. For the strains of Microcystis being investigated in the experiments, it seems increased algal growth with increasing temperature is responsible for the increase in colony size. In deriving our model, we have previously assumed aggregation due to cell growth is negligible, but this may not be true during peak surface water temperature conditions, leading to an underestimation of average colony diameter during high-temperature conditions. To account for cell growth in future iterations of this model, the method of Ackleh and Miller (2018) for calculating cell growth within a colony-where only a certain proportion of cells along the edge of the colony are able to reproduce new cellsshould be incorporated into Equation (7).
If we consider that quiescent conditions are hypothesized to be an immediate precursor to surface HABs (Michalak et al., 2013), then incorporating a growth term would likely change our results for intermittent high wind events (Figure 8c,d). We would expect slower frequencies of wind mixing to result in more opportunities for growth at the water surface during low wind conditions, leading to faster aggregation, which would cause a discrepancy between slow frequency and high-frequency wind mixing not currently demonstrated in this model. Recall that Yang et al. (2020) determined that intermittent disturbance not only promoted aggregation in M. aeruginosa but total biomass as well.  (Michalak et al., 2013). While this may be true immediately preceding bloom formation, it is also true that there must be enough large wind events before the quiescent period to encourage aggregation in order for a surface bloom to form. But, the occurrence of large wind events is still not enough: These wind events must occur when algal concentrations exceed 10 8 cells/m 3 in order for large colonies to form within a day. In addition to modeling concerns, this finding has implications for water quality management. If water samples are taken from well above the photic depth in a lake dominated by motile and colonial cyanobacteria, algal concentrations will likely be low and the average colony size will likely be quite small, which may give the appearance that HAB formation is unlikely. Meanwhile, large colonies could be rapidly forming at subsurface algal concentration peaks near the photic depth, indicating a surface bloom is imminent.

| Future work
A major objective of a mechanistic model is to generate hypotheses that drive further research. The results of this model suggest the need for a subsequent field study where meteorological conditions, lake thermal profiles, and both Microcystis concentration and colony size are tracked over depth and time at a relatively high frequency.
Once model results can be validated with field data, there are many further avenues of the study suggested by the model, both from an ecological and numerical perspective. One major ecological concern of M. aeruginosa is the ability to produce and release microcystins, a cyanotoxin. Microcystins are known to increase in extracellular concentration when Microcystis is stressed, and they also seem to have a relationship with extracellular polysaccharide content and colony size (Hu & Rzymski, 2019;Li et al., 2020;Rzymski et al., 2020;You, 2020). In fact, it is even hypothesized that microcystins can trigger colony formation via quorum-sensing processes (Rzymski et al., 2020). This raises an important question: How might the coupling of microcystin-triggered quorum sensing with colony dynamics improve model predictions of both the spatial heterogeneity of M. aeruginosa biomass and extracellular microcystin concentrations? After all, M. aeruginosa is a threat to public health because they release microcystins. In this regard, the fundamental question is not necessarily where the Microcystis is, but where the microcystins are.
Keeping in mind that the goal is to improve predictions over a seasonal time scale, then it will be necessary to use our model as a subroutine-in addition to a subroutine for disaggregation-in larger modeling software that can handle hydrodynamics, biogeochemical cycling, and algal life cycles (e.g., AEM3D (Hodges & Dallimore, 2016) or Delft3D-WAQ (Q. Chen & Mynett, 2006)). Since this model demonstrates aggregation is negligible except during high wind events at high algal concentrations, future models could also include a term that switches aggregation off when those conditions are not met. It would also be worthwhile to use these results to instead explore the evolution of the average colony size, d k , as a function of algal cell concentration and turbulence intensity. The model proposed in this manuscript is necessary to gain biological and physical insight into algal aggregation processes, but it may be possible to reduce some complexity once the system is understood.
Aggregation processes mostly affect buoyant transport, which is governed by the colony diameter-dependent settling velocity described in Equation (13). By restructuring the modeling in this way, the system of k equations can be avoided and bulk parameters remain the focus, removing most of the numerical expense that would be added by incorporating Equation (7) as a subroutine in software like AEM3D.
While the model described here has been derived for M. aeruginosa specifically due to their ubiquity and ecological importance, the modeling framework can easily be applied to any motile and colonial phytoplankton species. Different species have different motility and sticking mechanisms, so calculations of the advective velocity, w k (z,t), and sticking probability, α(i,j), will need to be tailored to each individual species. M. aeruginosa uses intracellular gas vesicles and buoyancy regulation mechanisms to achieve vertical motility, but many species of green algae use flagella to move about the water column, as an example. Despite these differences in subroutine calculations, the theoretical framework will remain largely unchanged from species to species and lake to lake. To promote the use of this model for different algal species, editable and annotated Matlab code used to simulate the base case scenario in Section 3.1 can be found at the Data Repository for the University of Minnesota (DRUM).

| CON CLUS ION
We have developed a theoretical model that tracks the meteorological-driven movement and aggregation of M. aeruginosa in lake water columns. There are limitations in this model-in particular, disaggregation is not accounted for and no explicit validation with field data has been made. However, the process of construct- Model results also clearly demonstrate that wind-induced mixing and algal concentrations exceeding 10 7 cells/m 3 are necessary to promote the aggregation of an initial single cell population to an algal population dominated by large colonies (d 101 = 160 μm) within 6 weeks, a time scale in accordance with field measurements (Cao & Yang, 2010;Xiao et al., 2018). This finding suggests quiescent conditions alone are not sufficient for surface bloom formation of colonial and motile harmful algae-large wind events prior to quiescent conditions are an important necessary precursor. In addition, the model provides guidance for future field data collection and model studies (e.g., quantifying the roles of extracellular polysaccharide and microcystin content as they relate to aggregate sticking probability). To practically implement the results of this theoretical model, we have identified ways to (i) incorporate this model into larger software in computationally efficient ways, and (ii) extrapolate this theoretical framework to different algal species.

ACK N OWLED G M ENTS
The authors wish to thank the Ramsey Lake neighborhood association, the St. Anthony Falls Laboratory engineering and tech support team, Dr. Shahram Missaghi, and Dr. Jiaqi You for help with the Ramsey Lake data acquisition.

CO N FLI C T O F I NTE R E S T
The authors cannot identify any potential conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data archiving is underway at the Data Repository for the University of Minnesota (DRUM), where interested parties can find lake thermal and hydrodynamic profiles from Ramsey Lake, MN, and an example Matlab simulation script.

O PEN R E S E A RCH BA D G E S
This article has earned Open Data and Open Materials badges. The data will be uploaded to the Data Repository for the University of Minnesota upon article acceptance. In the meantime, data has been included in the submission packet. Matlab code to run the model will be uploaded to the Data Repository for the University of Minnesota upon article acceptance. In the meantime, the code has been included in the submission packet. For all parameter values, see Table A1.
We assume that the cell density of each individual cell, whether it belongs to a colony or not, will react to the instantaneous light intensity in every grid cell at every time step by these equations.
To calculate colony density, ρ k , from single cell density, we use the relationship where n cell (k) is the ratio of cell volume to colony volume, n gas is the ratio of gas vesicle volume to colony volume, and ρ muc is the density of cell mucilage (Medrano et al., 2013). Like Medrano et al. (2013), we have kept n gas and ρ muc constant for all simulations; unlike Medrano et al. (2013), n cell is not constant but will vary with colony size because we assume all M. aeruginosa colonies have a fractal geometry and the relative content of mucilage to cells will increase with increasing cell size (Zhu et al., 2014). For parameter values, please see Table A1. Once ρ k has been calculated for all values of k at every time step in every grid cell, we can then calculate the buoyant velocity, w k (z,t), by Equation (13).